source('00_util_scripts/mod_bplot.R')
source('00_util_scripts/mod_bulk.R')
library(GEOquery)

proj.nm <- 'mission/SLE_TRPM2_MfMo/'

illu.m2 <- 'ILMN_2352380'

# IBD GSE169568 ---------
ibd1 <- getGEO('GSE169568', AnnotGPL = T) |>
  pluck(1)

ibd1@featureData@data |> as_tibble() |>
  filter(ID == 'ILMN_2352380') |> DT::datatable()

se1 <- ibd1 |>
  makeSummarizedExperimentFromExpressionSet(probeRangeMapper)

mean.exp <- tb1 |> filter(diagnosis.ch1 == 'Healthy control') |>
  summarise(mex = mean(exprs), .by = .feature)

mean.exp |>
  mutate(delta = abs(mex - 6.08955166666667)) |>
  slice_min(delta, n = 10)

tb1 <- ibd1 |> tidybulk()

tb1$diagnosis.ch1 |> table()

tb1 <- tb1 |>
  select(.feature, .sample, exprs, diagnosis.ch1)

res1 <- tb1 |>
  mutate(group = case_when(str_detect(diagnosis.ch1, 'Symp') ~ 'Symp',
                           str_detect(diagnosis.ch1, 'Heal') ~ 'HC',
                           .default = 'IBD')) |>
  filter(group != 'Symp') |>
  identify_abundant(factor_of_interest = group, minimum_counts = 5) |>
  test_differential_abundance(~ 0+group, contrasts = 'groupIBD-groupHC',
                              method = 'limma_voom',
                              omit_contrast_in_colnames = T) |>
  pivot_transcript()

m2.res1 <- res1 |> filter(.feature == 'ILMN_2352380') |>
  select(`logFC`, `adj.P.Val`) |>
  mutate(dataset = 'IBD GSE169568')

ibd1@featureData@data |> as_tibble() |>
  mutate(gene = ifelse(ID == 'ILMN_2352380', 'TRPM2', `Gene symbol`),
         .feature = ID, .keep = 'none') |>
  right_join(res1) |>
  select(gene, .feature, logFC, adj.P.Val) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/GSE169568.IBD.deg.csv')

m2.res1 <-
  read_csv('mission/SLE_TRPM2_MfMo/results/GSE169568.IBD.deg.csv') |>
  filter(gene == 'TRPM2') |>
  mutate(dataset = 'IBD_GSE169568')

read_csv('mission/SLE_TRPM2_MfMo/results/GSE169568.IBD.deg.csv') |>
  filter(gene != '') |>
  select(-.feature) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/IBD_GSE169568.deg.csv')

# IBD bulk GSE186507 -----------
ibd2 <- download_ncbi_counts('GSE186507')

ibd2.meta <- download_ncbi_meta('GSE186507', simple = F)

ibd2.meta |> dplyr::count(`diseasebi:ch1`)

tb2 <- ibd2 |>
  pivot_longer(-1)

tb2

tb2 <- ibd2.meta |>
  mutate(name = geo_accession, group = `diseasebi:ch1`, .keep = 'none') |>
  right_join(tb2) |>
  tidybulk(.sample = name, .transcript = Symbol, .abundance = value) |>
  quick_process_bulk(group = group)

tb2.abn <- tb2 |>
  filter(.abundant) |>
  select(-merged_transcripts)

res2.abn <- tb2.abn |>
  test_differential_abundance(~ 0+group, contrasts = 'groupIBD-groupControl',
                              method = 'edgeR_likelihood_ratio',
                              omit_contrast_in_colnames = T) |>
  pivot_transcript()

res2.abn |>
  filter(Symbol == 'TRPM2')

reseq2 <- tb2.abn |>
  test_differential_abundance(~ 0+group, method = 'DESeq2',
                              contrasts = list(c('group', 'IBD', 'Control')),
                              omit_contrast_in_colnames = T) |>
  pivot_transcript()

m2.res2 <- reseq2 |>
  filter(Symbol == 'TRPM2') |>
  mutate(dataset = 'IBD GSE186507', logFC = log2FoldChange, adj.P.Val = padj,
         .keep = 'none')

reseq2 |>
  select(Symbol, log2FoldChange, padj) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/GSE186507.IBD.deg.csv')

m2.res2 <-
  read_csv('mission/SLE_TRPM2_MfMo/results/GSE186507.IBD.deg.csv') |>
  filter(Symbol == 'TRPM2') |>
  mutate(dataset = 'IBD_GSE186507', logFC = log2FoldChange, adj.P.Val = padj,
         .keep = 'none')

read_csv('mission/SLE_TRPM2_MfMo/results/GSE186507.IBD.deg.csv') |>
  filter(Symbol != '') |>
  mutate(gene = Symbol, logFC = log2FoldChange, adj.P.Val = padj,
         .keep = 'none') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/IBD_GSE186507.deg.csv')

# psoriasis array GSE13355 ----------
pso1 <- getGEO('GSE13355', AnnotGPL = T)

pso1.m2p <- pso1[[1]]@featureData@data |>
  as_tibble() |>
  filter(str_detect(`Gene symbol`, 'TRPM2$')) |>
  pull(ID)

pso1.se <- pso1[[1]] |> makeSummarizedExperimentFromExpressionSet()

pso1.se |> select(contains('ch1')) 

pso1.se |> dplyr::count(characteristics_ch1)

pso1.tb <- pso1.se |> select(.feature, .sample, exprs, title) |>
  mutate_samples(group = str_extract(title, 'PP|NN|PN')) |>
  tidybulk()

pso1.reskn <- pso1.tb |>
  mutate(group = ifelse(group == 'NN', 'HC', 'PSO')) |>
  test_differential_abundance(~ 0+group, contrasts = 'groupPSO-groupHC',
                              method = 'limma_voom',
                              omit_contrast_in_colnames = T) |>
  pivot_transcript()

pso1.reskn <- read_delim('mission/SLE_TRPM2_MfMo/GSE13355.top.table.tsv')

pso1.reskn |>
  filter(Gene.symbol == 'TRPM2')

pso1.res <- read_delim('mission/SLE_TRPM2_MfMo/GSE13355.top.table (1).tsv')

pso1.res |>
  filter(Gene.symbol == 'TRPM2')

m2.res3 <- pso1.reskn |>
  filter(Gene.symbol == 'TRPM2') |>
  mutate(dataset = 'Psoriasis GSE13355')

pso1.reskn |>
  select(Gene.symbol, adj.P.Val, logFC) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/GSE13355.psoriasis.deg.csv')

m2.res3 <-
  read_csv('mission/SLE_TRPM2_MfMo/results/GSE13355.psoriasis.deg.csv') |>
  filter(Gene.symbol == 'TRPM2') |>
  mutate(dataset = 'Psoriasis_GSE13355')

read_csv('mission/SLE_TRPM2_MfMo/results/GSE13355.psoriasis.deg.csv') |>
  filter(Gene.symbol != '') |>
  dplyr::rename('gene' = 'Gene.symbol') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/Psoriasis_GSE13355.deg.csv')

# psoriasis array GSE109248 ------
pso2 <- getGEO('GSE109248', AnnotGPL = T)

pso2.se <- pso2[[1]] |> makeSummarizedExperimentFromExpressionSet()

## 17 Pso vs 14 HC
pso2.se |> distinct(title) |> DT::datatable()

pso2.tb <- pso2.se |>
  filter(str_detect(title, 'Control|Psoriasis')) |>
  select(.feature, .sample, exprs, title) |>
  tidybulk()

pso2.qtl <- pso2.tb$exprs |> quantile(c(0,.05,1))

pso2.res <- pso2.tb |>
  mutate(group = ifelse(str_detect(title, 'Pso'), 'PSO', 'HC')) |>
  keep_abundant(factor_of_interest = group, minimum_counts = pso2.qtl[2]) |>
  test_differential_abundance(~ 0+group, contrasts = 'groupPSO-groupHC',
                              method = 'limma_voom',
                              omit_contrast_in_colnames = T) |>
  pivot_transcript()

m2.res4 <- pso2.res |>
  filter(.feature == illu.m2) |>
  mutate(dataset = 'psoriasis GSE109248')

ibd1@featureData@data |> as_tibble() |>
  mutate(gene = ifelse(ID == 'ILMN_2352380', 'TRPM2', `Gene symbol`),
         .feature = ID, .keep = 'none') |>
  right_join(pso2.res) |>
  select(gene, .feature, logFC, adj.P.Val) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/GSE109248.psoriasis.deg.csv')

m2.res4 <-
  read_csv('mission/SLE_TRPM2_MfMo/results/GSE109248.psoriasis.deg.csv') |>
  filter(gene == 'TRPM2') |>
  mutate(dataset = 'Psoriasis_GSE109248')

read_csv('mission/SLE_TRPM2_MfMo/results/GSE109248.psoriasis.deg.csv') |>
  filter(gene != '') |>
  select(-.feature) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/Psoriasis_GSE109248.deg.csv')

# psoriasis array GSE41664 --------
pso3 <- getGEO('GSE41664', AnnotGPL = T)

pso3.m2 <- pso3[[1]]@featureData@data |> as_tibble() |>
  filter(`Gene symbol` == 'TRPM2') |>
  pull(ID)

pso3.se <- pso3[[1]] |> makeSummarizedExperimentFromExpressionSet()

## 24 pso vs 24 hc
pso3.tb <- pso3.se |>
  filter(characteristics_ch1.1 == '') |>
  mutate(group = ifelse(str_detect(tissue.type.ch1, 'Non'), 'HC', 'PSO')) |>
  select(.feature, .sample, exprs, group) |>
  tidybulk()

pso3.qtl <- pso3.tb$exprs |> quantile(c(0,.05,1), na.rm = T)
pso3.qtl

pso3.tb$exprs <- pso3.tb$exprs |>
  log2()

pso3.res <- pso3.tb |>
  mutate(exprs = ifelse(is.na(exprs) | exprs < 0, 0, exprs)) |>
  keep_abundant(factor_of_interest = group, minimum_counts = 0) |>
  test_differential_abundance(~ 0+group, contrasts = 'groupPSO-groupHC',
                              method = 'limma_voom',
                              omit_contrast_in_colnames = T) |>
  pivot_transcript()

m2.res5 <- pso3.res |>
  filter(.feature == pso3.m2) |>
  mutate(dataset = 'psoriasis GSE41664')

pso3[[1]]@featureData@data |>
  as_tibble() |>
  mutate(.feature = ID, gene = `Gene symbol`, .keep = 'none') |>
  right_join(pso3.res) |>
  select(gene, logFC, adj.P.Val) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/GSE41664.psoriasis.deg.csv')

pso3.res <- read_csv('mission/SLE_TRPM2_MfMo/results/GSE41664.psoriasis.deg.csv')

pso3.res |> filter(gene != '') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/Psoriasis_GSE41664.deg.csv')

m2.res5 <- pso3.res |>
  filter(gene == 'TRPM2') |>
  mutate(dataset = 'Psoriasis_GSE41664')

# IBD GSE119600 -------
ibd3 <- getGEO('GSE119600', AnnotGPL = T)

ibd3 <- ibd3[[1]]

ibd3@featureData@data |> as_tibble() |>
  filter(ID == illu.m2)

ibd3 <- ibd3[, str_detect(ibd3$`condition:ch1`, 'contr|Crohn|ulcera')]

ibd3$group <- ifelse(ibd3$`condition:ch1` == 'control', 'HC', 'IBD') |>
  fct_relevel('IBD')

ibd3.res <- geo_limma(ibd3)

ibd3.m2 <- ibd3.res |>
  filter(ID == illu.m2) |>
  mutate(dataset = 'IBD_GSE119600')

ibd3@featureData@data |>
  as_tibble() |>
  select(ID, `Gene symbol`) |>
  right_join(ibd3.res) |>
  dplyr::rename('gene' = 'Gene symbol') |>
  filter(gene != '') |>
  select(-ID) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/IBD_GSE119600.deg.csv')

ibd3.res |>
  mutate(logp = -log10(adj.P.Val)) |>
  tidyplot(logFC, logp) |>
  add_data_points()

# MS ----------
## GSE17048 ----
ms1 <- getGEO('GSE17048', AnnotGPL = T)

ms1 <- ms1[[1]]

ms1@featureData@data |>
  as_tibble() |>
  filter(ILMN_Gene == 'TRPM2') |>
  select(ID)

ms1$group <- ifelse(ms1$source_name_ch1 == 'blood, HC', 'HC', 'MS') |>
  fct_relevel('MS')

ms1.res <- geo_limma(ms1, use_vooma = F)

ms1.res |>
  write_csv('mission/SLE_TRPM2_MfMo/results/GSE17048.MS.deg.csv')

ms1.res |>
  filter(ID == 'ILMN_1743987')

ms1.se <- ms1 |> makeSummarizedExperimentFromExpressionSet()

ms1.se |>
  select(.feature, .sample, exprs, )

ms1.tb <- ms1.se |>
  select(.feature, .sample, exprs, group) |>
  as_tibble() |>
  mutate(exprs = ifelse(exprs < 0, 0, exprs) |> log1p()) |>
  tidybulk(.transcript = .feature, .sample = .sample, .abundance = exprs)

ms1.tres <- ms1.tb |>
  keep_abundant(factor_of_interest = group, minimum_counts = 0) |>
  test_differential_abundance(~ 0+group, contrasts = 'groupMS-groupHC',
                              method = 'limma_voom',
                              omit_contrast_in_colnames = T) |>
  pivot_transcript()

ms1.m2 <- ms1.tres |>
  filter(.feature == 'ILMN_1743987') |>
  mutate(dataset = 'MS_GSE17048', logFC, adj.P.Val, .keep = 'none')

ms1@featureData@data |>
  as_tibble() |>
  right_join(ms1.tres, join_by(ID == .feature)) |>
  select(`Gene symbol`, logFC, adj.P.Val) |>
  dplyr::rename('gene' = 'Gene symbol') |>
  filter(gene != '') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/MS_GSE17048.deg.csv')

## GSE41850 -----
ms2 <- getGEO('GSE41850', AnnotGPL = T) |> pluck(1)

fData(ms2) |> glimpse()

ms2.m2p <- fData(ms2) |>
  as_tibble() |>
  select(ID, gene_assignment) |>
  filter(str_detect(gene_assignment, 'TRPM2'))

# 23 HC vs 88 MS ?
pData(ms2) |> as_tibble() |>
  dplyr::count(`disease:ch1`, `treatment:ch1`, `visit:ch1`, `data set:ch1`) |>
  DT::datatable()

varLabels(ms2) <- varLabels(ms2) |> make.names()

ms2 <- ms2[, str_detect(ms2$treatment.ch1, 'untreat')]

ms2$group <- ifelse(ms2$disease.ch1 == 'Control', 'HC', 'MS') |>
  fct_relevel('MS')

ms2.res <- ms2 |> geo_limma(gene_col = 'gene_assignment')

ms2.res <- ms2.res |>
  separate_wider_delim(gene_assignment, '//', names = c(NA, 'symbol'),
                       too_many = 'drop', too_few = 'align_start') |>
  mutate(symbol = str_squish(symbol))

ms2.m2 <- ms2.res |>
  filter(symbol == 'TRPM2') |>
  mutate(dataset = 'MS_GSE41850')

ms2.res |>
  filter(str_starts(symbol, 'TRPM'))

ms2.res |>
  filter(symbol != '') |>
  dplyr::rename(gene = symbol) |>
  select(-ID) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/MS_GSE41850.deg.csv')

## GSE138064 -----
ms3 <- getGEO('GSE138064') |> pluck(1)

fData(ms3) |> glimpse()

fData(ms3) <- fData(ms3) |>
  separate_wider_delim(gene_assignment, '//', names = c(NA, 'symbol'),
                       too_many = 'drop', too_few = 'align_start') |>
  mutate(symbol = str_squish(symbol))

fData(ms3) |>
  as_tibble() |>
  filter(symbol == 'TRPM2')

# 8 HC vs 67 MS
pData(ms3) |> glimpse()

pData(ms3) |> select(title) |> DT::datatable()

ms3 <- ms3[, str_detect(ms3$title, 'Heal|0\\s{2}h')]

ms3$group <- ifelse(str_detect(ms3$title, 'Heal'), 'HC', 'MS') |>
  fct_relevel('MS')

ms3.res <- ms3 |> geo_limma(gene_col = 'symbol')

ms3.m2 <- ms3.res |>
  filter(symbol == 'TRPM2') |>
  mutate(dataset = 'MS_GSE138064')

ms3.res |>
  filter(symbol != '---') |>
  dplyr::rename(gene = symbol) |>
  select(-ID) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/MS_GSE138064.deg.csv')

# SS ---------
## GSE84844 -----
ss1 <- getGEO('GSE84844') |> pluck(1)

fData(ss1) |> glimpse()

fvarLabels(ss1) <- fvarLabels(ss1) |> make.names()

fData(ss1) |>
  filter(Gene.Symbol == 'TRPM2')

# 30 HC vs 30 SS
pData(ss1) |> glimpse()

varLabels(ss1) <- varLabels(ss1) |> make.names()

ss1$disease.ch1 |> table()

ss1$group <- ifelse(str_detect(ss1$disease.ch1, 'control'), 'HC', 'SS')

ss1.res <- geo_limma(ss1, gene_col = 'Gene.Symbol')

ss1.m2 <- ss1.res |> filter(Gene.Symbol == 'TRPM2') |>
  mutate(dataset = 'SS_GSE84844')

ss1.res |>
  filter(Gene.Symbol != '') |>
  dplyr::rename(gene = Gene.Symbol) |>
  select(-ID) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/SS_GSE84844.deg.csv')

## GSE48378 -----
ss2 <- pluck_geo('GSE48378')

fData(ss2) <- fData(ss2) |>
  separate_wider_delim(gene_assignment, '//', names = c(NA, 'symbol'),
                       too_many = 'drop', too_few = 'align_start') |>
  mutate(symbol = str_squish(symbol))

fData(ss2) |>
  filter(symbol == 'TRPM2')

pData(ss2) |> pull(disease.status.ch1) |> table()

# 16 HC vs 11 SS
ss2$group <- ifelse(ss2$disease.status.ch1 == 'healthy', 'HC', 'SS') |>
  fct_relevel('SS')

ss2.res <- geo_limma(ss2, gene_col = 'symbol')

ss2.m2 <- ss2.res |>
  filter(symbol == 'TRPM2') |>
  mutate(dataset = 'SS_GSE48378')

ss2.res |>
  filter(symbol != '') |>
  dplyr::rename(gene = symbol) |>
  select(-ID) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/SS_GSE48378.deg.csv')

## GSE51092 ---------
ss3 <- pluck_geo('GSE51092')

fData(ss3) |> glimpse()

fData(ss3) |> filter(Symbol == 'TRPM2')

# 32 HC vs 190 SS
pData(ss3) |> glimpse()

ss3$disease.state.ch1 |> table()

ss3$group <- ifelse(ss3$disease.state.ch1 == 'none', 'HC', 'SS') |>
  fct_relevel('SS')

ss3.res <- ss3 |> geo_limma(gene_col = 'Symbol')

ss3.m2 <- ss3.res |>
  filter(Symbol == 'TRPM2') |>
  mutate(dataset = 'SS_GSE51092')

ss3.res |>
  filter(Symbol != '') |>
  dplyr::rename(gene = Symbol) |>
  select(-ID) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/SS_GSE51092.deg.csv')

# RA ----------
disease <- 'RA'

## GSE15573 ------
acc <- 'GSE15573'

ra1 <- pluck_geo(acc)

fData(ra1) |> glimpse()

fData(ra1) |> filter(Symbol == 'TRPM2')

# 15 HC vs 18 RA
pData(ra1) |> glimpse()

ra1$status.ch1 |> table()

ra1$group <- ifelse(ra1$status.ch1 == 'Control', 'HC', 'RA') |>
  fct_relevel('RA')

ra1.res <- geo_limma(ra1, 'Symbol')

ra1.m2 <- ra1.res |> filter(Symbol == 'TRPM2') |>
  mutate(dataset = str_c(disease, acc, sep = '_'))

ra1.res |>
  filter(Symbol != '') |>
  dplyr::rename(gene = Symbol) |>
  select(-ID) |>
  write_source_csv(str_glue('{disease}_{acc}.deg'))

## GSE93272 ------
acc <- 'GSE93272'

ra2 <- pluck_geo(acc)

fData(ra2) |> glimpse()

fData(ra2) |> filter(Gene.Symbol == 'TRPM2')

pData(ra2) |> glimpse()

# 43 HC vs 232 RA
ra2$disease.state.ch1 |> table()

ra2$group <- ifelse(ra2$disease.state.ch1 == 'RA', 'RA', 'HC') |>
  fct_relevel('RA')

ra2.res <- geo_limma(ra2, 'Gene.Symbol')

ra2.m2 <- ra2.res |> filter(Gene.Symbol == 'TRPM2') |>
  mutate(dataset = str_c(disease, '_', acc))

ra2.m2

ra2.res |>
  filter(Gene.Symbol != '') |>
  dplyr::rename(gene = Gene.Symbol) |>
  select(-ID) |>
  write_source_csv(str_c(ra2.m2$dataset, '.deg'))

## GSE17755 ------
acc <- 'GSE17755'

ra3 <- pluck_geo(acc)

fData(ra3) |> glimpse()

fData(ra3) |> filter(Symbol == 'TRPM2')

pData(ra3) |> glimpse()

# 45 HC vs 112 RA
ra3$disease.ch1 |> table()

ra3 <- ra3[, str_detect(ra3$disease.ch1, 'indivi|^rheu')]

ra3$group <- ifelse(str_detect(ra3$disease.ch1, 'rheu'), 'RA', 'HC') |>
  fct_relevel('RA')

ra3.res <- ra3 |> geo_limma('Symbol')

ra3.m2 <- ra3.res |> filter(Symbol == 'TRPM2') |>
  mutate(dataset = str_c(disease, '_', acc))

ra3.res |>
  filter(Symbol != '') |>
  dplyr::rename(gene = Symbol) |>
  select(-ID) |>
  write_source_csv(str_c(ra3.m2$dataset, '.deg'))

# SLE -------
disease <- 'SLE'

## GSE154851 ------
acc <- 'GSE154851'

sle1 <- pluck_geo(acc)

fData(sle1) |> glimpse()

fData(sle1) |> filter(GENE_SYMBOL == 'TRPM2')

pData(sle1) |> glimpse()

# 32 HC vs 38 SLE
sle1$title

sle1$group <- ifelse(str_detect(sle1$title, 'Cont'), 'HC', 'SLE') |>
  fct_relevel('SLE')

sle1.res <- sle1 |> geo_limma('GENE_SYMBOL', force_normalize = T)

sle1.m2 <- sle1.res |> filter(GENE_SYMBOL == 'TRPM2') |>
  mutate(dataset = str_c(disease, '_', acc))

sle1.res |>
  filter(gene != '') |>
  select(-ID) |>
  write_source_csv(str_c(sle1.m2$dataset, '.deg'))

## GSE138458 ------
acc <- 'GSE138458'

sle2 <- pluck_geo(acc)

fData(sle2) |> glimpse()

fData(sle2) |> filter(Symbol == 'TRPM2')

# 23 HC vs 307 SLE
pData(sle2) |> glimpse()

sle2$case.control.ch1 |> table()

sle2$group <- ifelse(str_detect(sle2$case.control.ch1, 'SLE'), 'SLE', 'HC') |>
  fct_relevel('SLE')

# vooma not usable in non-complete cases
# but lmfit is ok
sle2.res <- sle2 |> geo_limma('Symbol', use_vooma = F, complete_case = F)

sle2.m2 <- sle2.res |>
  filter(Symbol == 'TRPM2') |>
  slice_min(adj.P.Val) |>
  mutate(dataset = str_c(disease, '_', acc))

sle2.res |>
  filter(Symbol != '') |>
  dplyr::rename(gene = Symbol) |>
  select(-ID) |>
  write_source_csv(str_c(sle2.m2$dataset, '.deg'))

## GSE65391 ------
acc <- 'GSE65391'

sle3 <- pluck_geo(acc)

fData(sle3) |> glimpse()

fData(sle3) |> filter(Symbol == 'TRPM2')

pData(sle3) |> glimpse()

# 23 HC vs 208 SLE
sle3$disease.state.ch1 |> table()
sle3$arthritis.ch1 |> table()
sle3$race.ch1 |> table()

sle3.aa <- sle3[, sle3$race.ch1 == 'AA']
sle3.aa$disease.state.ch1 |> table()

sle3$group <- ifelse(sle3$disease.state.ch1 == 'SLE', 'SLE', 'HC') |>
  fct_relevel('SLE')

# force_normalize greatly imporve significance in some cases
sle3.res <- sle3 |> geo_limma('Symbol', force_normalize = T)

sle3.m2 <- sle3.res |> filter(Symbol == 'TRPM2') |>
  slice_min(adj.P.Val) |>
  mutate(dataset = str_c(disease, '_', acc))

sle3.res |>
  filter(Symbol != '') |>
  dplyr::rename(gene = Symbol) |>
  select(-ID) |>
  write_source_csv(str_c(sle3.m2$dataset, '.deg'))

# sum-up --------
list(m2.res1, m2.res2, m2.res3, m2.res4, m2.res5) |>
  list_rbind() |>
  select(adj.P.Val, dataset, logFC) |>
  mutate(
    absFC = abs(logFC),
    direction = absFC/logFC,
    directed.p = -log(adj.P.Val) * direction) |>
  ggplot(aes(directed.p, dataset, fill = absFC)) +
  geom_col() +
  geom_vline(xintercept = c(-log(.05), log(.05)), linetype = 'dashed', color = 'grey') +
  geom_vline(xintercept = 0) +
  theme_pubr(legend = 'right') +
  scale_fill_gradient(low = 'blue4', high = 'red') +
  labs(title = 'TRPM2 expression from different autoimmune diseases',
       x = '-logPvalue', fill = 'log2FC') +
  theme_jpub

publish_source_plot('psoriasis.IBD.TRPM2.logfc', width = 70)

sum5 <- read_csv('mission/SLE_TRPM2_MfMo/results/psoriasis.IBD.TRPM2.logfc.csv')

list(sum5, ibd3.m2, ms1.m2, ms2.m2, ms3.m2, ss1.m2, ss2.m2, ss3.m2,
  sle1.m2, sle2.m2, sle3.m2, ra1.m2, ra2.m2, ra3.m2) |>
  map(\(x)select(x, adj.P.Val, dataset, logFC)) |>
  list_rbind() |>
  mutate(disease = str_extract(dataset, '^[:alpha:]+') |>
           fct_relevel('IBD','SLE','Psoriasis','RA','SS','MS')) |>
  mutate(disease_order = as.numeric(disease) * 1e6 + adj.P.Val, .by = disease) |>
  mutate(dataset = fct_reorder(dataset, disease_order, .desc = T)) |>
  mutate(
    absFC = abs(logFC),
    direction = absFC/logFC,
    directed.p = -log(adj.P.Val) * direction) |>
  ggplot(aes(directed.p, dataset, fill = absFC)) +
  geom_col() +
  geom_vline(xintercept = c(-log(.05), log(.05)), linetype = 'dashed', color = 'grey') +
  geom_vline(xintercept = 0) +
  theme_pubr(legend = 'right') +
  scale_fill_gradient(low = 'blue4', high = 'red') +
  labs(title = 'TRPM2 expression from different autoimmune diseases',
       x = '-logPvalue', fill = 'log2FC') +
  theme_jpub
  
publish_source_plot('AID6.TRPM2.logfc', width = 70)

## add covid & Flu ---------
aid_m2_logfc <- read_csv('mission/SLE_TRPM2_MfMo/results/AID6.COV.TRPM2.logfc.csv')

flu_m2_logfc <-
list.files('mission/SLE_TRPM2_MfMo/results/', pattern = 'flu',
           full.names = T) |>
  map(read_csv, id = 'file') |>
  purrr::reduce(full_join)

disea_m2_logfc <-
flu_m2_logfc |>
  mutate(adj.P.Val = ifelse(is.na(adj.P.Val), p_val_adj, adj.P.Val),
         logFC = ifelse(is.na(logFC), avg_log2FC, logFC),
         dataset = str_c('Flu_', (str_extract(file, 'GSE\\d+')))) |>
  full_join(aid_m2_logfc) |>
  full_join(j20_m2) |>
  mutate(disease = str_extract(dataset, '^[:alpha:]+') |>
           fct_relevel('IBD','COVID','SLE','Psoriasis','RA','SS','MS')) |>
  mutate(disease_order = as.numeric(disease) * 1e6 + adj.P.Val, .by = disease) |>
  mutate(dataset = fct_reorder(dataset, disease_order, .desc = T)) |>
  mutate(logp = -log(adj.P.Val))

disea_m2_logfc |>
  filter(!str_detect(dataset, 'CNP')) |>
  ggplot(aes(logp, dataset, fill = logFC)) +
  geom_col() +
  geom_vline(xintercept = -log(.05), linetype = 'dashed', color = 'grey') +
  geom_vline(xintercept = 0) +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  scale_fill_distiller(palette = 'RdYlBu',
                       values = pretty_distiller(disea_m2_logfc$logFC)) +
  labs(title = 'TRPM2 expression from different diseases',
       x = '-logPvalue', fill = 'log2FC') +
  theme_jpub

publish_source_plot('AID6.FLU.COV.TRPM2.logfc', width = 70)

disea_m2_logfc |>
  filter(!str_detect(dataset, 'CNP')) |>
  ggplot(aes(abs(logFC), dataset, fill = logp)) +
  geom_col() +
  #geom_vline(xintercept = -log(.05), linetype = 'dashed', color = 'grey') +
  geom_vline(xintercept = 0) +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  scale_fill_distiller(palette = 'Reds', direction = 1) +
  labs(title = 'TRPM2 expression from different diseases',
       fill = '-logPvalue', x = 'log2FC') +
  theme_jpub

publish_source_plot('AID6.FLU.COV.TRPM2.logfc2', width = 70)

# add sepsis ---------
aid6_m2 <-
read_csv('mission/SLE_TRPM2_MfMo/results/AID6.FLU.COV.TRPM2.logfc2.csv')

aid6_m2 <- aid6_m2 |>
  filter(!str_detect(dataset, 'RA|Psoria'))

sepsis_m2 <- list.files('mission/SLE_TRPM2_MfMo/results', full.names = T) |>
  str_subset('array.csv') |>
  map(read_csv) |>
  reduce(full_join) |>
  filter(gene == 'TRPM2') |>
  mutate(logp = -log10(adj.P.Val))

aid7_m2 <- aid6_m2 |>
  full_join(sepsis_m2)

aid7_m2 |>
  ggplot(aes(logp, dataset, fill = logFC)) +
  geom_col() +
  geom_vline(xintercept = -log10(.05), linetype = 'dashed', color = 'grey') +
  geom_vline(xintercept = 0) +
  scale_fill_distiller(palette = 'YlOrRd', direction = 1) +
  labs(title = 'TRPM2 expression from different diseases',
       x = '-log10(Pvalue)', fill = 'log2FC') +
  theme_jpub(theme_classic)

publish_source_plot('disease7.array.trpm2.logfc', width = 70)

aid6 <- list.files(pattern = '_GSE.+deg', recursive = T, full.names = T) |>
  map(read_csv, id = 'file')

tidy6 <- aid6 |>
  list_rbind() |>
  mutate(dataset = str_extract(file, '[:alpha:]+_GSE\\d+'), .keep = 'unused')

tidy6 |> dplyr::count(dataset)

# 376k <- 536k rows
tidy6 <- tidy6 |>
  slice_min(adj.P.Val, by = c(dataset, gene))

tidy6 |> write_source_csv('6aid.18datasets.deg')

kn.pal2 <- c('#d24040','#485ffc')

range6 <- tidy6 |>
  filter(str_detect(gene, '^TRPM[1-8]$')) |>
  pull(logFC) |> range()

pretty_distiller(range6)

tidy6 |>
  filter(str_detect(gene, '^TRPM[1-8]$')) |>
  ggplot(aes(dataset, gene, size = -log(adj.P.Val), fill = logFC)) +
  geom_point(shape = 21, stroke = .2) +
  scale_radius(range = c(0,3)) +
  theme_bw(base_size = 6) +
  theme(legend.key.size = unit(4, 'mm')) +
  rotate_x_text() +
  scale_fill_distiller(palette = 'RdBu',
                        values = pretty_distiller(range6)) +
  labs(color = 'log2FC', size = '-logPvalue',
       title = 'TRPM family gene expression from\ndifferent autoimmune diseases')

publish_source_plot('AID6.TRPM1-8.logfc.bubble', width = 70)

# HMCES ---------
aid6 <- list.files('mission/SLE_TRPM2_MfMo/results/', '_GSE', full.names = T) |>
  map(\(x)read_csv(x, id = 'file'))

aid6_hmces <- aid6 |> bind_rows() |>
  filter(gene == 'HMCES') |>
  mutate(dataset = str_extract(file, '[:alpha:]+_GSE\\d+')) |>
  slice_min(adj.P.Val, by = dataset) 

aid6_hmces |>
  ggplot(aes(-log10(adj.P.Val), dataset, fill = logFC)) +
  geom_col(color = 'black') +
  scale_fill_distiller(palette = 'RdBu',
                       values = pretty_distiller(aid6_hmces$logFC)) +
  theme_bw() +
  geom_vline(xintercept = -log10(.05), linetype = 'dashed') +
  labs(title = 'HMCES gene expression from different autoimmune diseases')
